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ABSTRACT 

Atmosphere is one of the most important noise sources for ground-based cosmic microwave back¬ 
ground (CMB) experiments. By increasing optical loading on the detectors, it amplifies their effective 
noise, while its fluctuations introduce spatial and temporal correlations between detected signals. We 
present a physically motivated 3d-model of the atmosphere total intensity emission in the millimeter 
and sub-millimeter wavelengths. We derive a new analytical estimate for the correlation between 
detectors time-ordered data as a function of the instrument and survey design, as well as several 
atmospheric parameters such as wind, relative humidity, temperature and turbulence characteristics. 

Using an original numerical computation, we examine the effect of each physical parameter on the 
correlations in the time series of a given experiment. We then use a parametric-likelihood approach 
to validate the modeling and estimate atmosphere parameters from the polarbear-i project first 
season data set. We derive a new 1.0% upper limit on the linear polarization fraction of atmospheric 
emission. We also compare our results to previous studies and weather station measurements. The 
proposed model can be used for realistic simulations of future ground-based GMB observations. 


1. INTRODUGTION 

Observations of the cosmic microwave background 
(GMB) are a unique probe of the fundamental physics at 
work in the early universe. Ongoing and upcoming high 
sensitivity observations of GMB temperature and polar¬ 
ization anisotropies provide constraints on the properties 
of the very early Universe, for example through observa¬ 
tions of GMB primordial polarization from large to inter¬ 
mediate angular scales. In addition, through the precise 
characterization of the impact of gravitational leasing on 
smaller scales, the analysis of GMB polarization brings 
information about large scale structures, and correspond¬ 
ing constraints on ne utrinos masses and species and on 
the dark energy, e.g., Planck Gollaboration et al. (2013); 


Das et al. (|2013). 

Ground-based experimental efforts have been playing 
and wil continue to play a prominent role in the exploita¬ 


tion of this potential (e.g. iThe POLARBEAR Gollabo 


ration et al.| 

2014); 

20T1); 

Naess et al.fl 

ipnS); WetlT (2 


BICEP2 Gollabora t ion et aT, . 
Crites et al.| ( |20I^ ; Abazajian et al 


instruments sucti as WlViAF[J or PlanclR or planned 
next-^neration space mission such as L)OrE[^ Lite- 
BIRE^and PIXIE ( |Kogut et ah] ( |2Qll[ )), ground-based 
experiments can deploy larger primary reflectors and 


planned 


^ http://map.gsfc.nasa.gov/ 

^ http://www.esa.int/Our.Activities/ 
Space_Science/Planck 

^ http://www.core-mission.org/ 

^ http://htebird.jp/eng/ 


therefore reach higher angular resolution, but must deal 
with residual atmospheric effects which can be an im¬ 
portant limitation to their ultimate performance. Emis¬ 
sion from water vapor and dioxygen molecules dominates 
atmospheric emission at millimeter and sub-millimeter 
wavelengths. This effect can be significantly minimized 
by observing in so-called atmospheric windows, where 
emission levels are greatly reduced. However, even this 
residual emission inevitably increases the optical power 
incident on the detectors and therefore their (photon) 
noise level, e.g., Arnold (2010). In addition, the inhomo¬ 
geneous distribution of water vapor molecules, driven by 
complex mechanisms, which depend on the properties of 
the atmosphere above a given observation site, results in 
both temporal and spatial variations of the received opti¬ 
cal power. The amplitude and correlation of atmospheric 
fluctuations depend on both the scanning strategy and 
the properties of the atmosphere at the time of obser¬ 
vation, such as the wind direction and speed. If treated 
as an additional noise-like component, this atmospheric 
contamination results in an additional colored and spa¬ 
tially correlated signal in the time stream of any specific 
detector. This would potentially affect not only the ef¬ 
fective noise level of the entire instrument, obtained as 
a combination of the noise levels of all its detectors, but 
also the rate at which this noise level decreases when the 
number of deployed detectors grows. 


polarized, e.g., Hanany & Rosenkranz 

( 

12003); 

Spinelli 

et al. (2011) and therefore appropriate 

hardware solu- 
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tions, such as the dual polarized TES (Arnold (2010)), 
may be used to mitigate its impact on the polarization- 
sensitive observations of the CMB. Nevertheless, even if 
such measures are implemented, the atmospheric emis¬ 
sion will in general contribute to the detected polarized 
signals due to instrumental limitations such as instru¬ 
mental polarization, imperfect h alf-wave plates or fre ¬ 
quency bandp ass mismatch (e.g., jShimon et al.| (|2008|); 
Errard| (|2012|)). Moreover, total intensity sensitivity 
IS crucial for calibration, which can impact, for exam¬ 
ple, the performance of comp onent separation based 
on a parametric approach, e.g. 

Therefore, though the focus of i. 


Stompor et al. (2009). 
ris work is on total in- 
tensity measurements, its conclusions are relevant for 
polarization-sensitive observations alike. 

Modeling of the atmospheric effects is complex. 
Eluctuations of the atmospheric optical depth generate 
emission with amplitude and scale which both depend 
on the observation site (dryness, air density, etc.) and 
on the time of observation (temperature, pressure, 
etc.). Moreover, wind can displace atmospheric struc¬ 
tures, introducing hard-to-mod el non-stationary effects . 
Nonetheless, it ha s been shown, Lay & Halverson (2000); 
Bussmann et al. (2005), that in the case of a given 
experimental setup (local plane, scanning strategy, 
etc.), and with help of some simplifying assumptions 
useful models can be developed and implemented. 
These previous efforts have used the 2d frozen screen 
approximation to characterize observations above the 
south pole. In contrast, we employ a 3d model, de rived 


as a generalization of the model of Church 


(1995) and. 


for the first time, apply it to the data gathered at the 
Atacama Desert. The modeling proposed here requires 
significant computational resources, but is more general 
and permits, for instance, accounting for near- and 
far-held beam regimes. 

The paper is organized as follows. In section we 
describe the total emission of the atmosphere, its im¬ 
pact on the photon noise of a ground-based detector as a 
function of frequency, and the typical emission law of 
atmospheric emission huctuations due to water vapor 
inhomogeneity. In section we derive a new general 
expression for the auto- and cross-correlation induced 
by atmosphere between two detectors of a given focal 
plane geometry. We present in section several results 
obtained with the model, computed using Quasi Monte 
Carlo algorithms, to illustrate the impact of the various 
parameters on the properties of the atmospheric signal. 
In section we compare the modeling predictions with 
real CMB data sets from the hrst season of polarbear- 
I observations. In addition, we derive estimates of the 
atmosphere model parameters, and compare them with 
measurements taken by a nearby weather station. We 
also compare our re sults with measurements made in Lay| 
& Halverson (2000) above the Atacama desert, and esti¬ 
mate an upper limit for the polarization fraction of the 
atmosphere. Einally, we conclude and discuss our re¬ 
sults in sectionAs a complement, we introduce in ap¬ 
pendix a faster method to estimate an approximation 
of the correlated noise between detectors, encapsulated 
in a full binned noise covariance matrix, and to efficiently 
simulate realistic noise time streams. 


2. ATMOSPHERIC TRANSMISSION AND 
EMISSION 

The atmosphere is only partly transparent at millime¬ 
ter wavelengths. Water vapor is responsible for most 
of the continuum absorption in the 100 GHz-1 THz fre¬ 
quency range. In addition, very strong absorption occurs 
at frequencies within the broad oxygen absorption band 
around 60 ± 10 GHz, as well as around the oxygen line at 
119 GHz and the water vapor lines at 22, 183, 325, and 
380 GHz. This leaves only four main atmospheric win¬ 
dows available for CMB observations from the ground: 
below 50 GHz, and around 95, 150 and 250 GHz. Ob¬ 
servations are also possible at about 340 GHz, but only 
when the atmosphere is exceptionally dry, i.e., for a low 
Precipitable Water Vapor (PWV) value. At higher fre¬ 
quencies, atmospheric opacity is incompatible with sensi¬ 
tive ground-based CMB observations, even from the best 
observing sites like the Atacama desert. 


PWV [mm] 

0.25 0.5 1.0 2.0 

time below 

5% 25% 50% 65% 


TABLE 1 

Percentage of time below the listed amount of total 

PRECIPITABLE WATER VAPOR, FOR THE LlANO DE ChAJNANTOR IN 

Chile. 


Table which summarizes the percentage of time at 
Atacama when the PWV is below a certain amount, 
shows that the typical amount of water vapor is about 
1 mm, the atmospheric conditions being better about 
50% of the tim^ We use the Atmospheric Transm is- 
sions at Microwaves (ATM) code (Pardo et al. (2001|)) to 
compute the atmospheric transmission 7 , detined as the 
frequency-dependent ratio I/Iq between the radiation / 
received by a detector and the radiation Iq above the at¬ 
mosphere. At the zenith, from the ALMA site located at 
5040 m elevation in the Atacama desert in Chile (Llano 
de Chajnantor), we get the transmission curves shown in 

Fig.[T] 



10 100 1000 
Frequency (GHz) 


Fig. 1. — Atmospheric transmission from the Atacama plateau at 
the zenith for different amounts of precipitable wate r vapor. This 
is obtained using the ATM code, [Pardo et aT] ( |200l[ ). 


Eor different zenith angles, the optical depth r has to 
be modified for the airmass m (90° — el), so that the at- 

^ https://almascience.nrao.edu/ documents-and-tools/ 
overview/ about- alma/atmosphere-model 
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mospheric transmission is 

T = exp(—r) = exp(—m (90° — el) tq). 


(1) 


The airmass m (90° — el) ranges from 1 at zenith (el = 
90°) to about 40 at the horizon (el = 0°), and tq is 
the optical depth at the zenith. The airmass can be 
computed for any zenith angle usin g the following fitting 
formula by Kasten & Young (1989): 


m{0) oc 


1 


cos 0 + 0.50572(96.07995 - 0 [deg])-i-6364' 


( 2 ) 


We note that for high elevations we have m(90° — el) 

1 / sin (el), as expected for a plane-parallel approxima¬ 
tion. 

The instantaneous noise of current CMB detectors is 
dominated by the statistical fluctuations in photons ar¬ 
riving at the detectors. Consequently instrument envi¬ 
ronment, specifically ground or space, plays a critical 
role. On the one hand, for space-borne detectors, the 
total background in the frequency range of interest for 
CMB observations is the sum of the CMB blackbody 
emission and the emission from the optics, dominated by 
the emission from telescope mirrors. 


I{iy) = Bjy{TcMB) + 


( 3 ) 


where Bj^(T) is the emission of a thermal blackbody at 
temperature T. For a space mission such as Planck, 
Ttei ^ 40 K, and on the basis of the pre-la unch measure¬ 


ments performed on the Planck reflectors (Tauber et al 
dM^ ), we can model the emissivity of a two-reflect or 


space-borne telescope as 


On the other hand, for ground-based experiments, the 
temperature of the telescope and cryostat window is of 
the order of 280 K and the total emissivity is typically of 
the order of a few percent. In addition, the emission of 
the atmosphere contributes to the total background with 
an additional term £{ 1 ^) of 

£{u) = [l-r{v)] 5,(Tatm), (5) 


where Tatm — 280 K is the atmosphere temperature. 
Thus, at 150 GHz, the total background emission is ^ 20 
K on the ground while it is ^ 1 K in space. 

In addition to this background emission, the emission 
of the atmosphere fluctuates as a function of time and 
pointing direction. A model of these fluctuations, due to 
inhomogeneities of its physical properties, is discussed in 
the next sections. 


3. ATMOSPHERE PHYSICAL MODEL 

In this section, we present the modeling of the atmo¬ 
spheric contamination occurring for ground-based ob¬ 
servations, through the derivation of an analytical ex¬ 
pression based on the auto- and cross-correlation be- 
tween detecto rs time streams, initially proposed by 


Church| (1995). Since the atmosphere mainly radiates 
non-polarized light, the detectors considered through¬ 
out the following sections are assumed to give measure¬ 
ments of the sky total intensity, i.e., the Stokes param¬ 
eter /, expressed in K. Eor example, in the case of 


antenna-coupled TES (Arnold (2010)), the considered 
time stream would effectively be the sum of the time- 
ordered data (TOD) coming from two orthogonal anten¬ 
nas confined within a focal plane pixel. 

3.1. Atmosphere eontrihution to antenna temperature 

Spatial inhomogeneities in the atmosphere emission 
can get imprinted in the TOD when the instrument line 
of sight scans across them. In addition, the wind blows 
structures through the beam. These inhomogeneities are 
due to turbulences along the line-of-sight, and are some¬ 
times assumed to be concentrated in a two-dimensional 
layer: atmospheric contamination is then approximated 
as a s creen of brightnes s temp e rature, moving with th e 


Sayers et al. 
are usually assume* 


wind (|Lay & Halverson! (12000 ) 

" (|201oj)-'^ - 

sumed t 


Bussmann et al. (|2005|) 
furthermore, these fluctuations 


to be frozen as the involved turbu¬ 
lent processes are mu ch smal l er tha n the displacements 
imposed by the wind (T aylor (|1938|)). 

The model proposed by Church 1(1995 ) is not restricted 
to a single turbulent layer, but instead sees the atmo¬ 
sphere as a continuum three-dimensional medium that 
depends on water vapor distribution, wind speed, tem¬ 
perature, etc. as one moves away from the telescope, 
along the line of sight. Our work consists in the imple¬ 
mentation and exploitation of this 3d modeling, hence 
significantly differing from previous s tudies. Based o n 
the Kolmogorov model of turbulence (|Tatarskii| (|1961|)), 
the power spectrum of the fluctuations in a large three- 
dimensional volume is 


P(k) oc k’’ 


(6) 


where k is the wavenumber with units and b = 

— 11/3. As a consequence, one can show that the con¬ 
taminated time streams have power law spectra, with 
an index relate d to the Kolmog o rov spectrum, typically 


P(f ) cc e.g., Marriage (2006) 


Errard 


(2012). Eollow- 


ing I Church (199’5), a coherence distance scale is defined 
over which tne Kolmogorov spectrum holds — the outer 
scale Lo corresponds to the typical distance that satisfies 
that condition. 

In the following reasoning, we use the geometry de¬ 
picted in Eig. An imaginary central detector of the 
focal plane observes along a vector fg = defined in 

spherical coordinates by (j) and 0 (equivalently azimuth 
and elevation). All other detectors are assumed to ob¬ 
serve in directions slightly off from the central line of 
sight, fg‘(t) = fs(^) + with 6P defined by the design 
of the focal plane. 

Each detector has an associated beam pointing in a 
given direction with an effective area H(fg^,r) such 
that, for a monochromatic detector observing at wave¬ 
length A: 


B{K,r) = 


where 


7rre^(f^ • r) 


X exp 


2(r^-(n-r)^) 

re^(f^ • r) 


-r) =wo\ 1 + 


X ri 


TTW^ 


( 7 ) 


(8) 


with wq is the beam waist given by 

A 


Wo = 


Trffb’ 


( 9 ) 
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Fig. 2.— A detector of the focal plane observes along a vector fg^ 
(line of sight), defined in spherical coordinates 0* and 6'^ (equiva¬ 
lently azimuth and elevation: 6^ = 0 corresponds to the horizon), 
fs is the pointing direction of an imaginary central detector onto 
the focal plane. The scan strategy, ss, can be expressed as the 
evolution of fg with respect to time, ss = drs/dt. Wind is a vector 
orthogonal to the 2 ;-axis and is characterized by its direction 
and its norm W. 



Fig. 3.— Illustration of three functions used in the modeli ng o f 
the correlation : the water vapor column, X 2 {z,z), Eq. JlSl; 
the temperature profile, Tphyg, Eq. ( |12| ; and the Gaussian beam 
width w{z)^ Eq. §. 


Oi) being the beam opening angle in radians. As expected, 
for large distances from the telescope, i.e., |fg • r| ^ rco, 
re(fg • r) has a slant asymptotic with slope wq, leading 
to 5(fg,r) oc 0^. It is computationally expensive to 
consider a realistic behavior for the beam effectiv e area. 


and this was usu a lly neglected in previous works (Lay fc 


Halver son (12000); Bussmann et ah (2005); Sayers et al 


(l-iOlOt). F ig7|3|shows the behavior of the Gaussian beam 
width w{z) given in Eq. compared to two other im¬ 
portant functions used in this modeling and introduced 
below, the water vapor density (Eq. (18)) and the tem¬ 
perature profile (Eq. ([T^). 

Given the expression for H(fg,r), the contribution 
dTant to the antenna temperature of a small element dV 
of atmosphere located at a distance r is proportional to 
the effective area of the telescope beam as seen from that 


point, i.e.. 


dT^ntir) = d a(r, A) Tphys(r) ^ (10) 

where the transmission T(A) mentioned in Eq. ^ is re¬ 
lated to the absorption a(r. A), expressed in m , by 


T{X) = 


(11) 


^phys(i*) in Eq. (10) is the physical temperature of the 
given volume of mmosphere. It is convenient to assume 
an adiabatic atmosphere, so that Tphys depends linearly 
on altitude z = r • z. 


^phys (^) ^phys ('2)) 


= T, 


ground 


1 - 




( 12 ) 


with Tground ~ 280 K the temperature at the ground level 
and 2;atm r\j 10^ m a typical height which depends on the 
observation site. Note that Eq. (12) is not expected to 


be correct above the tropo pause, generally located be 
tween 10 and 20 km, e.g. Hobgood (1993). We keep 


this assumption through the paper as our results do not 
depend significantly on 2(atm, the atmospheric emission 
being also weighted by the water vapor density profile. 
This latter decreases much more rapidly than Tphys with 
altitude as explained in the following paragraph. 


3.2. Analytical expression for the auto- and 
cross-correlation between detectors 


By integrating Eq. (10) over r, the total antenna tem¬ 


perature Tant as measured by a detector i can be written 
as 


^ant—^ant(rs (t)) 

1 r 

):2 ^(^s'(^)’b«(brphys(r). (13) 

More generally, using Eq. (f^, the correlation between 
two samples t and t', measui^ by two given detectors i 
and j, is defined as (Tg^nt( 0 ^ant(^ 0)5 expressed as: 


X ((a(r)(a(r )) Tphys (r) Tphys (r ) ( 14 ) 

where the average (•) is taken over realizations of the 
sky. As in|Church| (|1995|), we reduce the correlation (•) 
between two given points r and r' in the atmosphere to 
the ((a(r)a(r')) term. Note also that we neglected the 
wavelengths in ((a(r, A)a(r', A')). We will limit ourselves 
to the monochromatic case in this work but will mention 
the A 7 ^ A' possibility in section 

Quantities in Eq. ( |l4| ) depend on the atmosphere prop¬ 
erties (hidden in the (a(r)(a(r')) x Tphys(i*)Tphys(r') term) 
and on the experimental design and the operation of the 
telescope (included in the T(fg^(t), r)T(fgj(t'), r') term). 
At this point, if no wind is assumed, the only time de¬ 
pendence in Eq. ( [l4| ) is encoded in the scanning strategy 
i.e., fs(t). 
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= 100 m- 


To derive an analytic ex pression for th e correlation 
term, ((a(r)(a(r')), we follow Cliurcli| ( 1995|) and assume 
that 

(Q;(r)Q;(r')) = Xi{\r - r'\) X 2 {z, z') (15) 

where xi is the turbulent part of the correlation, effective 
for leng t hs sat isfying |r — r'| < Lo- In this work, as in 
Church (1995), xi is approximated by a Gaussian, i.e.. 


Xi(Ar) = x?exp - 


Ar^ 


(16) 


where Ar = |r — r'|. As depicted in Fig. this is nearly 
equivalent to Kolmogorov correlation, given by 

^Koim.(A^) cx: _ / {k Ar) hcdhz (17) 

Knn \ n 


where K^nin and /^max 


are proportional to the inverse of 
typical sizes of the turbulences, respectively their outer 
(energy provision) and inner (energy dissipation) scales. 
In the case depicted in Fig. the relative difference 
between the integration of Eqs. and 0 over Ar 
is about 10%. In the follow ing, we will use the Gaus¬ 
sian correlation given in Eq. (16) for computational effi¬ 


ciency. While the quantitative change due to the Gaus¬ 
sian approximation is less than 10% for , the Gaus¬ 
sian shape for xi drives the physical interpretation of the 
typical scale Lq. 


The second term in Eq. (15), X 2 ^ only depends on the 


altitude and can be interpreted as the water vapor dis¬ 
tribution. This latter is assumed to be an exponential 
decreasing with altitude z. 


X2{z,z')=X2exp ( - 


2 zq 


(18) 


which seems to be a reasonable approximation as in¬ 
dicated by some o bservations made abov e the Chilean 
desert of ^acama (Giovanelli et al.|(2001|)). Eor illustra¬ 
tion, Eig.[^ shows X 2 {z^z) compared to the beam width 
w{z) and me temperature profile Tphys(^)- The model 


used for X 2 naturally depends on the site of observa¬ 
tion. Eor example, water vapor density column at South 
Pole falls off qu ite significantly at altitu des outside a 300- 
2, 000 m layer ([Bussmann et al.| (|2005|)). We come back 
to this difference in appendix]A 

We assume in this work that wind behaves such that 
it does not mix atmosphere between different altitudes. 
Wind is thought to displace atmospheric structures, in 
particular the fluctuations encoded in the ((a(r)a(r')) 
correlation term. We suppose that wind affects the 
Xi(Ar) term of the correlation, Eq. (16), and leads to 
a new time dependence in this latter xi(Ar) -)> Xi(Ar') 
such as 


Ar' = Ar'{At) = Ar - WAt 


(19) 


with At = \t' — t\ and where W is a vector describing 
wind direction and amplitude. W is suppo sed to not 
depend on the a ltitude 2 ) as in Church (1995) or Lay & 
Halverson (2000), i.e.. 


W = lF^x + lF^y = lFw, 


( 20 ) 


where w corresponds to the unitary vector along the di- 
rection of the wind. It is possible to refine Eqs. ( [191^ 
and include a dependence of W on the altitude. We do 
not explore this possibility in this work and rather focus 
on a simple effective wind, described by a constant wind 
norm W = |W| [m s“^] and a single direction (pw [deg]. 
This approximation turns out to rea^nably describe real 
observations, as detailed in section 
Amplitude terms in Eqs. (16) anc 


dimension of C//' 

f'j 


IS 


ensure that the 
In this work, we set Xi = 1-0 


and X 2 — 1-0 The global amplitude of the corre¬ 

lation is parametrized by an effective ground tempera¬ 
ture, To, expressed in Note that measurements of 

the true ground temperature Tground, as measured by a 
weather station, and the estimation of Tq from , cf. 
section can give an estimate for the conversion factor 
between ATeff and K. 

3.3. Comparing the 3d and 2d approaches 


Previous studies such as Lay fc Halverson|(|2QQQ ); 


Buss¬ 


mann et al. (2005); bayers et al. T pTui r used th e approx¬ 
imation of a 2d frozen screen. Appendix [A] shows that 
the 3d modeling introduced above is general enough to 
reproduce the 2d screen predictions. In fact, parameters 
can be tuned so that the angular correlations from both 
approaches are nearly equivalent on small angular scales. 
However, this would assume a stepped function for the 
water vapor column, which is n ot supported by observ a- 
tions above the Atacama desert (Giovanelli et al.|(2001|)). 
Unlike measurements from south pole, estimated param¬ 
eters from a 2d screen approach might be difficult to in¬ 
terpret: in particular, the typical altitude h and thickness 
Ah of the turbulent layer would not have clear physical 
meanings in the case of distributed atmospheric turbu¬ 
lences. 

More generally, note that there are two significant rea¬ 
sons for using the 3d formalism. Eirst, the beams are 
overlapping where much of the emission is occurring, and 
this can be the case for large aperture instruments. Sec¬ 
ond, the scale length can be smaller than the thickness 
of the layer where the atmosphere is emitting. 
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Despite the caveats mentioned above, we will com¬ 
pare the effective brightness of the atmospheric ffuctu- 
ations, as estimated from both approaches. This com¬ 
parison is reasonable since the observed brightness of at¬ 
mospheric ffuctuations, should not depend on the 
assui ned model ing. At small angular scales, identifying 
Eqs. (Al) and (A6) leads in particular to Tq oc 5^, with 
a coefficient of proportionality depending on Lo, h. Ah, 
the angular scale 7 and ^^atm- In section we will have 
estimates for Tq and Lo derived from real observational 
data. Although significant uncertainties on h and Ah 
above the Atacama desert make any quantitative com¬ 
parison between the 2d and 3d approaches quite risky, 
estimates of atmosphere brightne ss d erived from both 
modelings are presented in section |5^ 


4. IMPLEMENTATION AND AN ILLUSTRATIVE 
EXAMPLE 

In this section, we illustrate the 3d model introduced 
in sectionby computing elements of Eq. ( [l4| , for 
various atmosphere conditions, as parametrized oy the 
atmospheric ffuctuations typical scale To, the wind W, 
the typical altitudes and Zg,tm^ as well as the effectiv e 
ground level temperature Tq. Contrary to Church (1995), 
no simplifications nor approximations are imposed on 
Eq. (14), and we perform the full computation of the 


six integrals to evaluate each C// element. 


4.1. Assumptions: scan strategy, focal plane layout 

Eor illustration purposes, we set a scanning strategy of 
the telescope, Ts{t), to be constant in elevation {Os{t) = 
Os^o) and sweeping in azimuth. As a rough approximation 
of usual scanning strategies adopted by CMB ground- 
based experiments, fs(^) is assumed in this section to 
follow a simple cosine function, 

cos {2tt 

fscant + 7/;), (21) 

where we set Aaz to be the angular size of the scan 
(usually of the order of a couple of degrees) and /scan the 
scan frequency, defined as 


/scan [H^-] 


ss [deg s 
Aaz [deg] 


( 22 ) 


where ss = |ss| is the norm of the scanning speed. Eor 
simplicity, we only consider a single detector in this sec¬ 
tion, pointing along the central line of sight rs(t) of the 
telescope. Throughout this paper, this detector will be 
denoted i = j = 0. Note that in section in which 
the modeled atmospheric contamination is compared to 
real data, Cft is computed using the real pointing of the 
experiment {(f = az{t),0 = el{t)}. 


change to the integral value. In addition, we use the spa¬ 
tial locality of the functions involved in the integrand to 
speed-up the computation, namely B{Ts{t),r), Eq. 
and Xi(Ar), Eq. ([T^ : we reduce the integration bounds 
in Eq. (14) to scales within the beam size or the turbu¬ 
lent typical length. Eor example, we reduce dr' to a 3d 
volume of radius 


fs(t) - r ± Lo, (23) 

which is effectively the la region of the correlations in¬ 
duced by the turbulences. In addition, the (6>, 0) and 
{0',(j)') variables are integrated around 


A max 



± max ( ^*6, L ) ) , (24) 


More distant locations are suppressed either by the beam 
term B{Ts{t), r), or by the turbulence term, xi- Outside 
the ‘’Icr” volume defined by Eqs. (23) and (24), inte¬ 
grand contributes to less than 5% oTthe totalintegral 


in C V . Einally, the implementation takes advantage of 


tt' 


symmetry relations such as C V = Cj- 
Eor a given set of physical parameters {Lo, Tq, ...}, the 
implementation of the algorithm used for this article per¬ 
forms the computation of a single element C// in ^ 1 


ms on one CPU. A basic parallelization of the algorithm 
is easy to implement and efficient to compute the entire 
set of C// elements. 

f'j 


4.3. Illustration of the modeling for a detector 
auto-correlation 

Assuming the specific observing strategy given in 
Eq. (21), Eig. [^depicts Cft' for the central detector 
i = j = 0 as a function of both t' {x-axis assuming t = 0) 
and one of the parameters involved in the atmosphere 
conditions description (^-axis, all the other parameters 
are set to their fiducial values). Cqo^ corresponds to 
the time-domain autocorrelation of the central detector. 
More precisely, Eig. ^depicts the quantity 

where is computed for the chosen fiducial pa¬ 

rameters { Lo = 300 m, zq = 2000 m, Zatm = 40,000 m, 
FWHM = 3.5arcmin, ss = 1 deg s~^, Aaz = 5 deg, 
050 = 0 deg, Os^o = 45 deg, -0 = 0 deg } and { IT = 25 
m s“^, <fw = 0 deg } in the windy cases (panels (5) and 
(6) only). A common feature across all panels in Eig.j^is 
the presence of nearly scan-synchronous signals, induced 
by atmospheric structures, roughly acting as a sky-like 
signal for low wind speed and/or short period of time. 
We describe below the effect of changing a parameter 
individually on the auto-correlation function: 


4.2. Implementation 

Significant computational resources are required to 
compute Eq. ( p^ . We use a Quasi Monte Carlo inte¬ 
gration methooto numerically estimate the six-integral 
included in the computation of a single Cft element. 
The number of quasi-random samples has been optimized 
and turned out to be ^ 10^ for the computation of a 
single C V element. More samplings lead to sub-percent 


• Turbulence typical length, Lq'. panel (1) shows 
that an increase of this length results in a larger 
width of the scan-synchronous features. In the 
limit of an infinitely large turbulent typical scale, 
all the detectors would be 100% correlated, for an 
infinitely long time if there is no wind. Contrarily, 
a Lo —> 0 m reduces the nearly scan-synchronous 
features, bringing the averaged correlation to neg¬ 
ligible levels. 
















7 




( 2 ) 

=300 m, 2 o =2000 m, Aaz =5 deg, 
W=0 m.s~^, <f>yy =0 deg 



( 5 ) 

=300 m, Zq =2000 m, ss = l deg.s"^, 
Aa2=5 deg, 0.17 =0 deg 



time [sec] 


( 3 ) 


Lq =300 m, Zq =2000 m, ss = 1 deg.s ^, 
Vr=0 m.s ~^, 0^7 =0 deg 



=300 m, Zq =2000 m, ss =1 deg.s ^, 
Aaz=b deg, M^=25 m.s~^ 



time [sec] 


1.35 


1.20 


1.05 


0.90 


0.75 


0.60 


0.45 


0.30 


0.15 


Fig. 5. — Normalized auto-correlation of a detector as a function of time, i.e., Cqq^ as defined in Eq. jl4| >, for various 

atmosphere physical parameters. The fiducial parameters are chosen as { Lo = 300 m, VF = 0.00 m s“^, c^w = 0-00 deg, Aaz = 5.00 deg, 
ss = 1.00 deg s~^ }. The observing experiment is assumed to scan the sky in azimuth (here 6 back-and-forth movements in 30 s), at a 
given elevation. A detailed description of each panel can be found in the main text. 


• Correlation amplitude, Tq: as mentioned ear¬ 
lier, this is an effective ground temperature and its 
square is the global normalization of the correlation 
in the case of monochro- 
Tq simply scales the 


matrix , Eq. 
matic detectors, 
entire correlation function, and is not illustrated in 

Fig.|5j 


(|14), 

Eerefore 




Telescope scanning properties, ss and Aaz: 
they impact /scan, Eq. (22), and consequently the 
frequency of the scan synchronous features in . 


As shown in panel (2), in the absence of wind, the 
larger the scan speed is, the more scan-synchronous 
features are present. As depicted in panel (3), the 
larger Aaz is (for a fixed scan speed), the fewer 
scan-synchronous features appear. A long and slow 
(respectively short and rapid) movement of the 
telescope modulates the atmospheric structures to 
low (respectively high) time streams frequencies. 


• Typical altitudes z^tm and zq: they have simi¬ 
lar impact on the correlation. Panel (4) of Fig. 
shows the effect of changing zq on Cqo^ : smaller 
water vapor column results in smaller correlation 
amplitude, with a mild impact due to the form of 
X 2 , cf. Eq. (18). In addition, note that the position 
of the water vapor densest region would also play 
a role on the width of the nearly scan-synchronous 
features. The X 2 term acts as a weight for the 
contribution of turbulent structures to the atmo¬ 


spheric signal: the lower is zq, the more important 
is the contribution of atmospheric fluctuations at 
lower altitudes. 


• Wind properties: results are depicted on pan¬ 
els (5) and (6). The first one shows the effect of 
wind speed on Cqo^ , for a constant wind direc¬ 
tion (^w = 0 deg, cf. the geometry illustrated in 
Fig-E The larger is the wind speed, the more 
rapidly harmonics of Cqo^ are suppressed. Since 
wind is displacing atmospheric turbulences, any de¬ 
tector cannot observe the same sky features at a 
given azimuth pointing. From Eq. (14), having a 

non-zero W leads to C//' oc 


getting about the cross-term in this expression, this 
is why the correlation is suppressed for long time 
intervals \t' — t\ ^ LojW. 


Panel (6) shows the effect of a wind direction 0 vf 
change at a given wind speed, W = 2bms~^. 
The geometry of the scanning strategy with the 
wind direction is depicted in Fig. as defined 
in Eq. ( [^ , the assumed scanning strategy is cen¬ 
tered around {(j){t))t = 0 deg since we took the fidu¬ 
cial (j)s^o = 0 deg. On one hand, pseudo-harmonic 
structures of Cqq^ are efficiently suppressed for 
(j)w ^ { — 180, 0, +180} deg. For the chosen scan¬ 
ning amplitude Aaz = 5 deg, these wind directions 
are roughly orthogonal to the scanning direction, 
which means that the wind efficiently brings new. 






































Fig. 6. — Scan strategy of the telescope is a cosine function, 
Eq. ( |21| >, centered on (f)s,o = 0. The wind is assumed to be parallel 
to the X — y plane and have a fiducial direction (f)\Y = 0. 


and hence uncorrelated, atmospheric fluctuations 
across the telescope line of sight. On the other 
hand, even if atmospheric structures are shifted 
and suppressed as a function of At, significant 
features remain in the time streams for the case 
(pw = ±90 deg. The physical interpretation of 
this effect is that the telescope resamples already 
scanned structures on the sky: the pw = ±90 deg 
wind shifts atmosphere turbulences in a nearly par¬ 
allel direction to the telescope scans. Specifically, 
given our assumptions for the scanning strategy, a 
(pw = — 90 deg would move atmospheric structures 
with the line of sight motion during the 1st, 3rd, ..., 
(2i^^±l) subscan of the telescope. On the contrary, 
a pw = ±90 deg would move fluctuations with the 
line of si ght motion during the 2 nd, 4th, ..., 


subscan. 


Bussmann et al. 


description of this phenomenon 


(2005) gives a detailed 


The remarks above hold for cross-correlations, i.e., for 
i 7 ^ j case. However, we always have < CprJ and 

the wind can impact correlations amplitude and shape 
between detectors, depending on their specific positions 
on the focal plane, as well as on the scanning strategy. 

Complementary to Fig. which illustrates the cor¬ 
relation in the time domain. Fig. shows the Fourier 
transform of several Cqq , i.e., the auto-spectra, for 
various wind speed and wind directions. First, the upper 
panel of Fig. shows the effect of various turbulence 
scales Lo at a fixed wind direction and speed (VF = 25 
m s“^ and pw = 0 deg): similarly to the observations 
derived from Cqq , the auto spectra present features 
appearing at /scan = 0.2 Hz. Large Lq results in 
more power at low frequencies while small Lo leads 
to a broader spectrum, with spread power at high 
frequencies. Second, the middle panel shows the effect 
of various wind speeds at a fixed wind direction and Lo 
{pw = 0 deg and Lo = 300 m). The larger the wind 
speed, the broader are the peaks in frequency space. 
As mentioned previously, wind displaces atmospheric 
structures and the detector does not resample these 
structures at a constant frequency /scan- Third, the 
lower panel shows the effect of various wind directions at 
a fixed wind speed and Lo". as suggested by the correla¬ 
tion behavior in time domain, wind direction modulates 
the typical frequencies of atmospheric contamination 
features. 



0.0 0.1 0.2 0.3 0.4 0.5 


frequency [Hz] 

Fig. 7. — Auto power spectra, obtained as the Fourier transform 
of Cqq^ , shown here for various Lq (upper panel, assuming VF = 25 
m“^ and bvv = 0 deg), various wind speeds (middle panel, assum¬ 
ing 0VV = 0 deg and Lq = 300 m) and various wind directions 
(lower panel, assuming VF = 25 m s“^ and Lq = 300 m). 


To summarize, we presented the dependence of the 
autocorrelation regarding the turbulence typical size, on 
the observing strategy, on the wind properties and on the 
water vapor column typical height. In the limit of a low 
wind speed, the model predicts the presence of nearly 
scan synchronous features in the detector time streams. 
However, the width of these features strongly depends 
on the typical size Lo of the turbulence, and depends 
weakly on the typical altitudes zq and z^tm- Moreover, 
the amplitude of these features usually decreases rapidly 
as a function of the elapsed time — t|, mostly driven 
by the wind speed. The pseudo-periodicity of these 
features is modulated by the wind direction relative to 
the scanning strategy. 

The next section will quantitatively compare the 
predictions of the modeling described above with re¬ 
cent observations of the CMB ground-based experiment 
POLARBEAR-L 


5. QUANTITATIVE COMPARISON BETWEEN 
THE MODELING PREDICTIONS AND REAL 
CMB DATA SETS 

In this section, we fit real Polarbear-i data sets with 
the modeling introduced in section and derive distri¬ 
butions of atmosphere parameters such as typical turbu¬ 
lence scal e and wind speed . Our approach here is quite 
similar to Bussmann et al. (2005|), but differs in two as¬ 
pects: first, we use the 3d modeling of the atmosphere 
and second, we estimate the physical atmospheric pa¬ 
rameters through the maximization of a parametric like¬ 
lihood based on the correlation matrices C/- . A quanti- 
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— CES#i 


/ CES#I 


a new set of atmospheric parameters (except zq and 
^atm as explained below) to be estimated and therefore 
a new likelihood to be maximized. 

For a given polarbear-i CES, we estimate the full 
real datasets’ covariance matrix as 


el 

CES#0 


CES#N 


D 


tt' 




An,S) 

^3 


{m, n} G { t , t '} 


S 


(25) 


horizon, el=0 


az 



Fig. 8.— Drawing of the scanning strategy of ~ 8-hour 
POLARBEAR-l observations. It is composed of several constant ele¬ 
vation scans (CES). 


tative comparison with previous results derived fro m the 


2d modeling approa ch above the Atacama desert (Lay 
& Halverson (2000)) is also presented. This analysis 


involves an estimator of the full cov ariance of the real 
data sets introduced in section |5.1[ and a param etric 
maximum-likelihood framework detailed in sect ion 15.2[ 
Results of the fit are summarized in section lOl 


5.1. Analysis of POLARBEAR-I data sets 

POLARBEAR-l is a CMB polarization imaging experi¬ 
ment!^ located on the Chajnantor plateau, at 17, 000 feet, 
near to (and hence looking through a similar atmosphere 
as) the Atacama Cosmology Telescop^ the ALMA|^ 
and the APE}(|^ observatories. It started its operations 
in early 2012 and has performed observations of three 
patches of the sky at 150 GHz, leading to an effective 
coverage of 30 deg^ of the sky. The polarbear collab¬ 
oration has recently published the results of the first 
season of observation (2012-2013), and demonstrated 
measur ements of CMB H-modes at sub-degree angular 
scales ( Ade et al. (|2014]); IPOLARBEAR Collaboration 
et al.| (2013|); j'Lhe POLARBEAR Collaboration et al. 

p!l ) — - 


POLARBEAR-I is composed of 1274 bolometric detectors 
(637 dual-polarization pi xels sensitive to t o tal int ensity) 
sampled at 190.73 Hz (jKermish et al.| (|2012|)). As 
schematically depicted in h'ig. polarbear-i CMB 
observations are composed of constant elevation scans 
(denoted CES hereafter), each of them lasting ^ 15 
minutes. Each CES is itself composed of ^ 100 subscans 
i.e., right-to-left and left-to-right going scans, with a 
Acj) 5 deg amplitude. The motion of the telescope 
is a bit more complex than the one considered in 
Eq. (21), with consecutive acceleration/deceleration 
(turnarounds) and constant velocity periods. The scan¬ 
ning strategy (6>s(t), 0s (t)) is an input to our modeling, 
and we use the true motion of the Polarbear-i tele¬ 
scope in this section. The instrument periodically scans 
in azimuth and, as assumed in section for a given CES 
we are in a ^s = constant case. In this work, we fit each 
CES independently: each of them is considered as an 
independent realization of the atmospheric fluctuations. 


® http://bolo.berkeley.edu/polarbear/ 
http://www.princeton.edu/act/ 

® http://www.almaobservatory.org/ 

® http://www.apex-telescope.org/ 


where denotes the calibrated total intensity data 

from a pixel i of the focal plane, at a given time sam¬ 
ple m —sample which comes from a group of subscans 
S. The average {')s is taken over several periods made 
of consecutive subscans. In each of these groups, we 
consider an even number of subscans, so that left-going 
(right-going) subscans are always averaged with left¬ 
going (right-going) subscans only. 

Given the instrument specifications, the complete 
T)*f would have a size (upixeis x risampieses)^ ~ (10®)^. 
It would therefore be challenging to quickly perform 
operations on this data set. To simplify the computation 
of the real data covariance matrix D// , we use some 
properties of the atmospheric contamination. Eor 
reasonable scanning strategies, since the atmosphere is 
evolving at rather low temporal (< 1 Hz) and spatial 
frequencies (> 0.5deg), a few detectors among the focal 
plane with low frequency time streams can capture 
enough information about atmospheric fluctuations 
properties. We hence choose “cardinal” pixels across 
the focal plane (e.g., at the center and edges), and we 
downsample the time streams of each of these pixels in 
order to lighten DfJ and therefore speed up the analysis 
without degrading the parameters estimation. 


Eirst the compression of {m, n} into temporal bins 
{t, t'} in Eq. (25) corresponds to the downsampling pro¬ 
cess. Second, in Eq. (j^ is also the result of an 
average of the TOD over some regions of the focal plane, 
precisely over detectors close to the “cardinal” pixels i.e.. 




{dn 


k G neighbor pixels to pixel i ’ 


(26) 


where the neighbor pixels k are the nearest observing 
detectors to a “cardinal” pixel i. Complementary to 


Eq. (25), the error bar on each D// , denoted 


tt ' 


is taken as the standard deviation of 


\ ^ ! {m, n} G { t , t '} 

computed over groups of subscans S. We compute a 
single D matrix for each polarbear-i CES. 


Atmospheric contamination encapsulated in DP can 
sometimes be approximated by a stationary process, 
meaning that DP with At = t' — t. This 

pseudo-stationarity is exploited in Appendix where 
an approximation of is introduced for rapid charac¬ 
terization of colored and correlated noise. However, wind 
usually breaks stationarity since it modulates the atmo¬ 
spheric signal periodicity, which depends on the pointing 
of the telescope and on the direction of the wind. Eor 
example, in the case of a wind nearly parallel to the 
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scanning strategy {(pw = iOO deg in Fig. [^, the de¬ 
tectors of a left going focal plane or of a right going focal 
plane woul d not show the same correlation pattern, cf. 
section 14.31 


the physical parameters p = {Lo, Tq, W,...} describing 
the atmosphere conditions, with the real-data matrices, 
, estimated from the polarbear-i first season data 

sets. 


5.2. Parametric likelihood and implementation 

We want to construct a likelihood to quantitatively 
and robustly compare the predictions of the atmospheric 
contamination modeling, QP = CP (p), function of 


The comparison of modeled and real covariances has 
been tackled in other astrophysical contexts, for example, 
in the case of astrophysical foregrounds separation from 
multi-frequency data set s. Following the formalism by 


Pham & Cardoso (2001), the negative log-likelihood C 
can be written as 


- 2 log (>C(p)) oc 




log 


det C 


,*'(p) (d//) 


-1 


— constant > 


E 1*4(c “'(p) - D «') (D,f )■' (cy (p) - D «') (Dy)”'] |, 


(27) 

(28) 


where Eg. (p8|) is the quad ratic approximation of Eq. (27) 
(Delabrouille et al. ( 2003[ )). For the computation of 
note that QP (p) is given by the integration of Eq. (14) 
for a given set of parameters p, cf. sectionand DP is 
the covariance matri x com puted from polarbear-i real 
data sets, cf. section [5T| 

Eor the estimation of D, we consider time streams 
d-^ over groups of 6 consecutive subscans, defining S 
in Eqs. (25) and (26): averages {')s are then taken over 
all the a\imable groups of 6 consecutive subscans within 
a CES, corresponding to ^ 35 s long time streams. As 
mentioned earlier, note that S is composed of an even 
number of subscans, so that left -going scans are only av - 
eraged with left-going scans, cf. Bussmann et al. (2005). 
Similarly, right-going scans are only averaged over right- 
going scans. Prior to Eqs. (25) and (26), we remove the 
global mean signal from all the TOD over the course of a 
given CES. A ground-template, obtained as the estima- 
tio n of an azimuth-fixed signal, is removed as descr ibed 
by iThe POLARBEAR Collaboration et al.| (2014|). In 
addition, we choose 5 cardinal pixels (one at the cen¬ 
ter and four on the edge, forming a square shape). The 
average in Eq. (26) is taken over the 5 closest working 
pixels to each oithe cardinal pixels. We downsample 
the time streams from 190.73 Hz to 2 Hz, in order to 
speed up the algorithms without missing the main fea¬ 
tures of the atmospheric conta minat ion. As shown in 
Eig. and described in section |5.3| these features are 
essentially pseudo scan-synchronous occurring at a fre¬ 
quency of ^ 0.25 Hz in the case of polarbear-i. 

In addition, we provide the optimization algorithm 
with the derivatives of the likelihood, with respect to 
the considered parameters, namely djCjdLo^ dCjdW^ 
dC/d(pw and dCjdTQ. These derivatives are semi- 
analytically computed using Eq. ( [2^ , and include deriva¬ 
tives of the modeled covariance CP defined in Eq. (14). 
Having semi-analytical derivatives of the likelihood leads 
to a significant speed-up of the optimization. 


We implement and test the algorithms on NERSC 
systemf^ and use a simple parallelization scheme in 
which a single process performs the analysis of one 
CES; this corresponds to the computation of a DP and 
the estimation of the parameters {To, W, (pw, Tq} by 
maximizing the likelihood given in Eq. ( [2^ . Eor this 
work, we decide to fix zq = 2000 m and Zatm = 40, 000 
m as they turn out to be badly conditioned: the de¬ 
pendence of GP on these quantities is quite weak, and 
make the atmospheric model dege nerate, in particular 
with the amplitude Tq (section |4.3[ ) . We therefore allow 
only {To, VF, pw^ To} fo vary. In addition, scan speed, 
elevation and azimuth involved in Eq. ( p!4| ) are set by the 
real observations of the telescope, i.e., we compute the 
modeled C for the detectors location used to compute 
D, using the same scanning strategy, and at the same 
binned samples 


Einally, as illustrated in Eig. the estimation of the 
atmospheric modeling parameters is performed in two 
steps: we first compare with precomputed tem¬ 

plates G^P^, estimated for various atmospheric condi¬ 
tions and saved on disk. The minimization of —21og(>C) 
gives a rough estimate of the adjusted parameters, which 
are used as a starting point for the second step. We maxi¬ 
mize the likelihood expressed in Eq. (28), based this time 
on the matrices DP and GP , using gradient informa- 


tion in a trun cated Newton algorithm, e.g., [Wright fc 


Nocedal (2006). 


5.3. Results 

Eig. shows the correlations between the five 
“cardinal” pixels, + constant ± (colored 

points), and the corresponding estimated predictions 
from the atmospheric model, G^^P^ + constant (solid 
lines), with constant = —min{D^^P^). Only four 

National Energy Research Scientific Computing Center ma¬ 
chines http://www.nersc.gov/ 
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Fig. 9.— Schematic illustration of the algorithm we used to estimate the parameters of the atmospheric modeling from the ra w d ata. 
We first filter a ground template from the raw total intensity time streams, and then compute the correlation following Eq. J25l. A 


second step consists in comparing with precomputed (lighter on disk), leading to a rough estimate of the optimal parameters. 

This is the starting point for a minimization of the likelihood given in Eq. ( |28| >, over constrained variables p = {To, wind speed, wind 
direction. To}, using gradient information in a truncated Newton algorithm, c.f. Nash (1984), and based on the full matrices and 

Finally, we compute the final from the estimated atmospheric parameters. As an illustration, comparison of several and 

their adjusted are shown in Fig, 




ig’l^ 


parameters are adjusted here, namely the turbulence 
typical size (I/o), the wind speed and direction (W) and 
the effective ground temperature (Tq). The four panels 
have been picked from among the analysis of the entire 
first season of polarbear-i, corresp onding to 8432 CES 


for the three polarbear-i p atches (The POLARBEAR 
Collaboration et ah ( 2Q14[ )). The depicted CEb are 


taken on different days at different elevations (cf. title on 
each panel), and are chosen to show various atmospheric 
conditions. 


As predicted by the model, nearly scan-synchronous 
features are present in the real data sets, with a de¬ 
creasing amplitude as a function of time. As mentioned 
in section the amplitude of the curves depends on 
the effective temperature Tq, the slope of the envelope 
is a function of the wind speed, the width of the peaks 
is related to Lq and their relative positions depends 
on the wind direction. These results show that a 
four-parameter physical model of the atmospheric 
contamination seems able to describe the full noise 
covariance matrix estimated from real CMB data sets. 

A set of {I/O, VE, 0VF: To} is estimated for each CES. 
Eig. pr] shows the distributions of each of the estimated 
parameters, for all the analyzed CES. Eor this particular 
POLARBEAR-I data set, Lo, wind speed, wind direction 
and temperature distributions seem to have broad peaks 
around ^ 300 ± 100 m, ^ 25 ± 10 m s“^, ^ Odeg and 
^ 60 ± 30 mICeff respectively. The wind direction turns 
out to be weakly constrained from the POLARBEAR-I 
data sets: as mentioned in the previous section, this 
parameter only affects the relative position of the 
scan-synchronous peaks, and those can be quite small 
in some cases. Moreover, the small amplitude Aaz of 
the subscans does not offer a strong lever arm on (pw 
The peaks around (j)w = 0 deg and (j)w = 360 deg 
in the wind direction histogram correspond to the 


lower and upper bounds imposed during the likelihood- 
optimization, indicating that the algorithm is unable to 
find better optimal coordinates given the provided data 
sets. 

A comparison between Tq as estimated from the CMB 
data set and the true ground temperature as measured 
by the POLARBEAR weather station gives a rough 
estimate for the conversion factor between ICeff and K. 
An average of 0°C over the year at the ground level 
leads to the relation 1 mKQf^ ^ K. 

An interesting quantity derived above is the typical 
size of turbulence, Lq. The distribution of estimated 
values shows a peak around 300 m, which is larger 
than the assumed 1 — 100 m scale in Church (1995) 
and lower than the observe d 10 km ab ove the Owens 
Valley Radio Observatory in Lay ( 1997|). However, the 
estimation of Lq depends strongly on tne site location, 
and measurements seem to vary significantly between 
diffe r ent fr equency range of observation. In particular. 


Lay (1997) points out that interferometers on Mauna 
Kea or Plateau de Bure in Erance have observed quick 
phase variations, indicating the presence of rather small 
fiuctnation scales. Realistic atmosphere simulations 
in the optical reg i me ab ove the Atacama desert, e.g., 


Masciadri et al.| (2013), indicate that atmospheric 
ffuctuations can take place down to the 100 — 500 m 
scales, which is in agreement with our estimate of Lq. 


To evaluate the goodness of fit of the likelihood- 
optimized models presented in Eig. we compute the 
following 


- Cil’ip) 


= 




ad; 


tt' ( 


(29) 


where p corresponds to the likelihood-optimized atmo¬ 
sphere parameters. We show in Eig. the distribu- 
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11/25/2012, <el> = 43.26deg 
L^=129m, cp^=42.6deg 
W=4.09m/s, To=212.50mK^,, 


12/25/2012, <el> = 70.44deg 
L^=317m, cp^=0.15deg 
W=13.6m/s, To=70.25mK^^ 
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Fig. 10.— Example of estimated elements (colored points with error bars given by defined in Eq. ( |25| . They are shown 

here as a function of At = \t' — t\, each color corresponding to a different pixels pair {i, j}. Solid lines are the best fit model 
obtained by optimizing Eq. ( |28| over four atmospheric parameters: the turbulence typical length (To), the wind speed and direction (W 
and (pw) cind the effective ground temperature (Tq). These solid lines tend to overlap in some cases, but are always in good agreement 
with the observational data sets. Each panel corresponds to an independent CES, and the best fit parameters are detailed in the top-right 
corner of each figure, together with the date of the observation and the average elevation. 


tion of the reduced computed for each CES using 
Eq. (29). This distribution peaks around 1 as one could 
expect for models sensibly describing the real observa¬ 
tions. We should point out that the four free parameters 
{I/O, W, 0 VF 5 ^ 0 } might be slightly degenerate between 
each other. Nevertheless, error bars on the estimation of 
these parameters, computed usin g th e second derivative 
of the likelihood at the peak, Eq. (28), give uncertainties 
of the order of 5-10% of the central values. Comple¬ 
mentary to the we also evaluated the probability to 
exceed (PTE). We find that 25% (8%) of the CES have 
a reduced with a PTE lower than 5% (1%). Eig. 10 
shows cases which have a reduced < 3.0. The signifi¬ 


cant part of the data which has a low PTE can be caused 
by atmospheric conditions that break our assumptions— 
in particular the hypothesis of frozen turbulences or the 
assumption that wind has a flat profile with altitude. 


5.4. Comparison with Lay & Halverson (200C^ results 
above the Atacama desert and with weather 
station measurements 


We quantitatively compare in this paragraph the pre¬ 
dictions of our new modeling with independent measure¬ 
ments taken by a weather station, and with previous 
works which assumed the two-dimensional frozen screen 
approach for the atmospheric emission. 

Eirst, as mentioned in appendix we can compare the 
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Fig. 11. — Distribution of several physical parameters, recovered using the procedure described in section]^ Lq (left panel), wind speed 
{ sec ond) , wind direction (third) and effective temperature (right panel). These histograms are derived from the likelihood optimization, 
Eq. (|28|>, over each CES of the first season of polarbear-i. 



Fig. 12. — Distribution of reduced over all the analyzed CES, 
evaluated using Eq. ( |29| >. The dashed vertical line corresponds to 
a reduced x^ of 1. 


estimation for t he brightness of fluctuations, 

using Eg. (|A 6 |) with the measurements from [Lay fc 


Halverson (jzDCFO). This latter presents measurements 
from an interferometer observing the sky at 11 GHz (con¬ 
verted in the article to a brightness at 40 GHz) in the 
context of the 2d frozen screen approximation. This in¬ 
terferometer is composed of two 1.8 m-dishes, and we will 
assume that the instrument observes 7 ^ 1 deg angular 
scales on the sky. Since this latter angle and polarbear 
probed angular scales are both rather small, a compari¬ 
son between Tq and is possible but has to be treated 
carefully, cf. section |3.3| and appendix The alti tude 
and thickness of the turbulent layer being uncertain, |Lay| 
& Halverson ( |2QQQ[ ) considers a range of possible values 


for B'f, corresponding to various values of h (altitude of 
the turbulent layer). Even though this latter work quotes 
amplitudes of the fluctuations in \mK‘^ m], w e 

follow here the notations from Bussmann et al. (2005), 


with a brightness taken as 




(30) 


which is valid in the limit of lon g averaging time. Table 
summarizes the values quoted in Lay & Halverson (2000 
extrapolated from 40 GHz to 150 GHz such as 


^150 GHz —^40 GHz 


HSOGHz 


, (31) 

P40GHZ y 

c. {Al.ltlli) SloGHz for PWV = 1 ± 1 mm, 

where r is the opacity introduced in Eq. Table also 
shows the inferred values for the fluctuations brightness 
from POLARBEAR- 1 data sets, using estimations of Tq and 
Eq. (A 6 ). Both sets of quoted results have significant 
error bars, associated with the conver sion from refrac¬ 
tive in dex to brightness temperature in |Lay fc Halverson] 
( 200 0 |), a nd due to the variations of Lq in our study, cf. 
Eq. (A 6 ). In addition to providing a new semi-analytical 


comparison between 2d and 3d atmosphere modeling, we 
show that both sets of measurements turn out to be la 
compatible. 

Second, to quantify the correlation between the atmo¬ 
spheric parameters estimated from the GMB data sets 
and the ones coming from the POLARBEAR site weather 
station measurements, we compute the Spearman rank- 
order correlation coefficient, p, between the two data sets. 
We only use the best constrained data sets, i.e., such that 
their reduced < 3.0, Eq. ( [^ and Eig.[^ In this way, 
we remove the worst PTE cases. We flnd^ = 46.8% for 
the wind speed, p = 35.5% for the wind direction and 
p = 65.2% for the ground temperature. Because the 
estimation gives the atmospheric properties at an alti¬ 
tude where the atmospheric fluctuations are located, the 
correlation between the weather station and the GMB 
data sets is not expected to be very strong. In fact, we 
find only a slight correlation between the two data sets. 
To complement this comparison, we show in Eig. p^the 
wind speed values as measured by the weather station 
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Lay & Halverson 

(2000 

1 

This work— 

from To, Lq 

and Eq. ( 

A6 

I 

Quartile 

500 m 

1000 m 

2000 m 

500 m 

1000 m 

2000 m 


25% 

50% 

75% 

1 oq + 128 

^-457 

238412489 

1808Rt|i 

78091*123 

1 900 + 1286 

59601«199 

2466012*6|o 

259 ± 89 

985 ± 337 

3377 ± 1155 

488 ± 179 

1854 ± 678 

6358 ± 2324 

820 ± 357 

3112 ±1355 

10671 ± 4647 


TABLE 2 

Brightness of th e atmospheric fluctuations, in mK‘^ above the Atacama desert, Chile. Comparison of results from|Lay & | 
I Halverson| (|2000|) (scaled from 40 CHz to 150 CHz using Eq. ([^) and this work, as a function of the co nsidered altitude 

FO R THE TURBULENT LAYER h, AND AS A FUNCTION OF THE QUARTILE OF MEASUREMENTS DISTRIBUTION. SIMILARLY TO |LAY & HALVERS0N| 

i 20001) , WE ASSUME A THIC KNESS A/l = h FOR THE TURBULENT LAYER. THERE ARE 50% UNCERTAINTIES ON THE QUOTED RESULTS FROM 
^AY & HALVERS0N| (|2000|), ASSOCIATED WITH THE CONVERSION FROM REFRACTIVE INDEX TO BRIGHTNESS TEMPERATURE. We COMBINE 
THIS ERROR WITH THE UNCERTAINTY ASSOCIATED TO THE CONVERSION FROM 40 TO 150 CHZ DUE TO THE PWV VALUE, EQ. ( [^ . We 
CONSIDER A TURBULENT SCALE Lq = 300 M FOR THE ESTIMATION OF THE BRIGHTNESS FROM OUR ESTIMATES OF Tq, AND INCLUDE ± 100 M 
VARIATIONS AROUND THIS VALUE (CF. LEFT PANEL OF FlC. LEADS TO ~ 30% UNCERTAINTIES ON OUR RESULTS. 



Fig. 13.— Comparison of the wind speed as measured by the weather station at the ground level (vertical axis) with the atmospheric 
wind as estimated from the polarbear-i first season data sets analysis (horizontal axis). Radii of circles are proportional to the inverse 
of the reduced Eq. ( |29| , and orange line corresponds to x = y. The side panels correspond to the normalized distributions for the 
measurements taken by iTTe weather station (right panel) and for the wind speed estimated by our analysis (top panel). 


at the ground level with the wind speed estimated from 
the POLARBEAR-l data sets analysis. Along with the dis¬ 
tributions of each data sets, data points are shown as 
circles, with sizes proportional to the inverse of their re¬ 
duced This illustrates similar distributions of wind 
speeds across the year, but with noticeable differences: 
(1) the wind speed as measured from is on aver¬ 

age higher than the ground wind speed and (2) the slope 
of a possible linear fit would not be equal to 1. Both 


these remarks might be explained by the fact that wind 
at the ground level is not necessarily the same as the 
one displacing the atmospheric turbulence. This latter 
may take place at a different altitude, with another am¬ 
plitude and potentially direction. A more relevant cor¬ 
relation would have to be performed between the atmo¬ 
spheric parameters estimated from CMB data sets and 
the measu rements of a balloon pr obe above the observa¬ 
tory, as in Bussmann et al. (2005). In addition to giving 
sensible priors for the likelihood optimization, Eq. (28), 
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altitude atmospheric characterization would provide cor¬ 
rect forms for the water vapor column density function, 
Eq. (18), the temperature, Eq. (12), and the wind profile. 


5.5. Upper limit on the polarization fraetion of 
atmosphere emission 


We perform the same analysis described in section 5.1 


but with data correlation matrices D/f estimated from 


pixel difference time streams instead of pixel sum. In 
such case, the nearly scan-synchronous features described 
earlier are highly suppressed and the global amplitude of 
the correlation is significantly lower. The absence of fea¬ 


tures in Tyfj makes the optimization of the likelihood 
more challenging, and usually results with a large wind 
speed and a low Tq. Eig. [14] shows the distribution of the 
amplitude Tq from the analysis of pixel sum time streams 
(i.e., total intensity data, similar to the right panel of 
Fig. [TT| and from the pixel difference time streams (i.e., 
polarized data). We compare these results with the es¬ 
timation of To from simulated white noise TOD, taken 
as random Gaussian realizations assuming a distribution 
width equals to the standard deviation of the real time 
streams. 

This simple approach shows that the amplitude of the po¬ 
larized atmospheric signal is at least two orders of magni¬ 
tude below the total intensity signal, from 60 ± 30 
down to < 300 /iT'eff • Moreover, we verify that this latter 
value is compatible with the analysis of simulated white 
noise time streams. This indicates that constraints on 
To from polarized data set mainly come from the uncor¬ 
related detector noise. Hence, from the analysis of the 
POLARBEAR-l first season of observation, the linear po¬ 
larization fraction of the atmosphere emission p satisfies 
p < 1.0%. 



Fig. 14.— Distribution of the fitted amplitude of the correla¬ 
tions, To, obtained from the analysis of t he t otal intensity (gray 
histogram, similar to the right panel of Fig. [nj, polarized data set 
(orange histogram) and from white noise simulated TODs (blue 
histogram). 


6. CONCLUSIONS 

We present a new modeling for the atmospheric emis¬ 
sion in the millimeter and sub-millimeter wavelengths. 


and fit it to the Polarbear-i data sets. This work 
introduces the equations to evaluate the six-dimensional 
spatial correlation induced by turbulent water vapor. 
This correlation is a key object to realistically simulate 
a three-dimensional atmospheric medium, and also 
simulate the TOD of a detector associated to a given 
scanning strategy. By comparing the model predictions 
with POLARBEAR-I data sets, we derived distributions 
of the main atmospheric parameters such as the typical 
size of the turbulence, wind speed, temperature, etc. 
The semi-analytical equations, along with the typical 
distributions of physical parameters, are the ingredients 
to perform realistic simulations of atmosphere. This 
ability might be crucial to evaluate and optimize the 
performances of future ground-based CMB experiments. 


We first review in section the atmospheric transmis¬ 
sion and background emission in the range of frequencies 
involved in CMB observations. We describe the thermal 
loading induced by atmosphere as well as the frequency 
and PWV dependence of the average contamination. 

In addition to the loading due to the atmosphere, the 
fluctuations present in this turbulent medium induce a 
correlated contamination for the detectors time streams 
of a given focal plane. On the grounds of the model¬ 
ing origmally introduced by Church (1995), we derive in 
section a new expression for both the auto- and cross¬ 
correlation induced by atmospheric fluctuations between 
two detectors of a given focal plane geometry. 

We present in section [^several results from an original 
numerical computation of the model, performed using 
Quasi Monte Carlo algorithms. The computation of the 
modeling involves a six-dimensional integral for a given 
pair of detectors and a given time sample: numerical ap¬ 
proximations of the integrals are required, and turn out 
to give sensible results. The main prediction of the mod¬ 
eling is the presence of nearly scan-synchronous features 
in the correlation function, modulated by the scanning 
strategy, the wind speed and direction, the atmospheric 
fluctuations typical size and the atmosphere tempera¬ 
ture. These features are however only significant for a 
limited set of atmosphere parameters, and may not be 
relevant for the estimation of astrophysical or cosmolog¬ 
ical quantities. 

We compare, in section the prediction of the mod¬ 
eling with real CMB data sets from the polarbear-i 
first season of observation, using a maximum para¬ 
metric likelihood approach. We derive estimates of 
four atmosphere parameters, {To, Tq, VE, (fw}^ and 
show that their distribution over the year were peaked 
around {300 ± 100 m, 60 ± 30 25 ± 10 m s“^, 

0 deg} respectively. The weakest parameter estimation 
concerns the wind direction (/)w, which mainly impacts 
the relative position in time of the features present in the 
cross-correlations. The most likely value, (j)w = 0 deg, is 
interpreted as an artifact of the analysis: it corresponds 
to a bound imposed on the parameter during the 
likelihood optimization. We do not see a significant 
improvement by cutting out cases with high reduced 
We compare our estimate of Tq with the brightness of 


atmos pheric fluctuations measured by [Lay fc Halverson 
(2000) above the Atacama desert. To perform this 
comparison, we use a semi-analytical 2d layer approxi¬ 
mation of the angular correlation estimated from the 3d 
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modeling; of the turbulences. Measurements by |Lay fc 


Halverson (2000) and from this work are la compatible 


k'inally, we give an updated observational upper limit 
on the polarization fraction of the atmospheric emission 
at 150 GHz: from the analysis of the polarbear-i 
polarized data sets, we find that the linear polarization 
fraction is < 1.0%. 


External information about wind, temperature and 
water vapor profiles as measured by e.g., a flying 
weather probe above the considered observatory would 
certainly help for the estimation of the atmospheric 
conditions from the CMB data sets. Measurements of 
turbulence characteristics performed by an independent 
instrument such as MASS/DIMM would certainly be 
useful too (Kornilov et al. ( 20Q7|) ). Adding priors to 


the likelihood equation, Eq. (|28|), can potentially speed 
up its maximization and more robustly constrain the 
different modeling parameters. In addition, note that 
we treat here each polarbear-i CES independently 
for the parameters estimation, but one could imagine 
grouping several CES together (for which atmosphere 
properties would be assumed to be stable), possibly 
helping in estimating or constraining parameters. 


Parametrization and estimation of can be of great 
potential for the analysis of ground-based experiment 
data sets. In particular, the proposed modeling gives 
a basis for the development of simulation tools. We pro¬ 
pose to study in a following work the effect of systematic 


effects on the performance of CMB polarization observa¬ 
tions, in the presence of realistic atmospheric contami¬ 
nation, detector noise and polarized sky signal. Coming 
CMB experiments, such as Stage-IV instruments, will use 
hundreds of thousands detectors: understanding, charac¬ 
terizing and potentially treating correlated noise will be 
critical for a recovery of their full sensitivity and an op¬ 
timized exploitation of their data sets. In particular, we 
note that having multichromatic detectors, and poten¬ 
tially atmosphere monitoring detectors outside the at¬ 
mospheric windows (although still in the optically thin 
limit), can help data analysts to better estimate the at¬ 
mospheric signal, using for example frameworks similar 
to t he astrophysical for egrounds separation techniques. 


Leach et al. (^2008). We also leave the implementa- 


e.g. _ _ 

tion and testing of this possibility to a future work. 
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CYT). Einally, we would like to acknowledge the tremen¬ 
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A: SEMI-ANALYTICAL COMPARISON WITH THE FROZEN 2 D SCREEN APPROXIMATION 


The approximation of a 2d frozen screen assumes an angular correlation given by (Lay & Halverson (2000)) 


C2d 


r 

,( 7 ) = 27r / 

Jo 


ip P{'ip) Jo (27r'07) 


(Al) 


where t/) is the angular wave number, P(V’) is the angular power spectrum and Jo is the Oth-order Bessel function. It 
is assumed that the observed field 7 is small and that the turbulent atmospheric layer is high i.e., 


Ah > 


h 


2 sin( 6 >) 


7- 


(A 2 ) 


This layer is located at an altitude h with a thick ness Ah, and the observ ation is made at an elevation 0. In the case 
of the Kolmogorov turbulence power law, Eq. ([^, [Bussmann et ah (2005) gives an expression for P{'ip) of the form 




(A3) 


where the small-angle approximation (7 <C 1 ) is employed. 

We show in the following derivation that the 3d modeling is general enough to recover the angular correlation of the 
2d screen approach. To compare Eq. (Al) with the correlation CpJ introduced in Eq. (14), we can collapse the 3d 

iphfication that the lines-of-sight for detectors i and j are pointing 
(1995), Eq. (14) can be expressed as 


modeling into a thick 2d layer. We assum e for sin 


at the zenith, i.e., rg, oc z. Similarly to Church 




C'mJoM ^ I dZx 2 {Z) {T^i.ys{Z)f X 


w'^jz, 7 ) 
2Ll 


-1 


(A4) 


where Z = {z z')/2^ and we assum ed that w‘^{z) 'Zi w‘^{z') which is only valid if \z — z'\ Z. We also use the 
notation X 2 {Z) = X 2 {z^z'). Eq. (A4) is derived by expressing Eq. (p^ in cartesian coordinates, and by changing the 


integration variables from {x, i, x , yP z'} to {(x -h x')l2^ (y + y')/z^ {z -\- z')/2^x — x'^y — yP z — z'}. The integral 
over Z can be limited to a turbulent 2d layer by taking a step function for X 2 such as 

Ahl 


X 2 {z,z') = 1.0m ^for{z,z'} G 


h± 


(A5) 
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-^h±^ 


We can hence derive an approximation of the collapsed 2d correlation model ( 7 ) from Eq. ( |A4| |, 

\2N -1 


C^3d^Idel(7)=^V^i^oT2(^l + 


{w{h, 7 ))^ 


X [1.0 m ] X A/i I 1 — 


h 


‘ ' V ^atm, 

where we took a constant Gaussian beam width w{Z) ~ w{h) ~ h^, and assumed that Ah <C .Zatm- 


(A6) 


Eq. (A 6 ) shows that larger Lq leads to larger correlation for a fixed angular scale 7 . A thicker {Ah /^) or a warmer 
(To Zg,tm /^) atmospheric layer would increase the correlation between measurements separated by an angle 7 . 
Moreover, larger angular separation 7 would lead to smaller correlation for fixed T^, To, h and Ah. Fig. [Tsl shows the 

h±^ — 

dependence of the correlation functions C2d screen ( 7 ) from Eq. (Al) and model (^) Eq. (A 6 ), as a function 

of the angular scale 7 . We look at various h G {500 m, 1000 m, 2O00m|, take = 100 m and Ati = 500 m. We 
consider a lower bound h/(2Ah) for the integral in Eq. (|A 1 |), since this sets the range of -0 for which Kolmogorov 

IJ ^^ h±^ 

(2000)). In all the considered cases, the fractional difference between model( 7 ) 

o deg. 


spectrum is valid (Lay & Halverson 
and C2d screen ( 7 ) is < 10 % tor 7 < 


This shows that the 3d modeling is general enough to reproduce the 2d screen correlation function. In fact, some 
parameters such as {h. Ah, T^, ^^atm} can be tuned so that C2d screen — m^odei small angular scales. 


B: APPROXIMATION OF THE FULL NOISE COVARIANCE MATRIX 
We propose an algorithm to build an approximation of DjU', less expensive computationally, lighter on disk, easier 
to handle. In the approximation that DjU is stationary, its Fourier transform is expected to be diagonal. We therefore 
build the full binned covariance matrix, named Dj^ in the following: i and j are individual detectors (or equivalently 
sum and differences of perpendicular detectors within a focal plane pixel) and 6 is a given frequency bin. However, 
note that the stationarity assumption is broken in reality and the formalism below could not be used to recover e.g., 
the wind properties from real data sets. 

We first introduce the algebra and its implementation, and second propose a framework to generate TOD from . 


Binned noise eovarianee estimation from TOD 

Similarly to what was proposed in Chiang et al. (2010), each considered time stream is apodized and Fourier 
transformed: 


'it e S, =FFT {Wt X d7) 


(Bl) 


where S denotes a group of consecutive subscans and Wt is the normalization of a hanning window, Wt, such that 


Wt = 


Wt 

(Wt%' 


(B2) 


Windowing prevents aliasing between frequency modes, which would bias the evaluation of the noise, especially at low 
frequencies. To describe the structure of the correlations between detectors as a function of frequency, we choose a set of 
^bins frequency bins 6 , which should be optimized to encapsulate noise properties (scan synchronous contaminations, 
low frequency power, high frequency effects, etc.). We describe below a suitable choice for this set, in the case of 
POLARBEAR. For a given CES, the full binned covariance matrix is defined as 


d4 = 




(f, s) 


feb 


(B3) 


Dk is a complex object of size n^et x '^det x Ubins: which satisfies the symmetry relation 

DA = (d7* (B4) 

In Fig. Il^ are depicted several typical shown here as 2-dimensional objects (in the {det i, det j} space), 
for three different bins b. In this example, i and j indices correspond to single detectors, i.e. polarbear-i 
bolometers. Because of the stationarity assumption, Dmisses some information about atmospheric contamination. 
However, one can observe the detector-detector correlations, which can be quite significant at low frequencies, / < I Hz. 


To avoid mode coupling and to better estimate the low frequency part of the noise, we consider groups of 25 
consecutive subscans for the analysis. Note that S can be composed of an odd number of subscans since the phase 
information (left- and right-going scans are not equivalent because of the wind) is anyway destroyed during the 
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Fig. 16.— Three normalized shown as 2d objects, for different frequency bins b. The x and y axis correspond to individual 
detector indices: in this case, i and j correspond to individual bolometers. The applied normalization is the following: \/i,j,b, —)► 

X The depicted quantity has therefore no units and the matrices diagonals are equal to 1. One can notice that the 

correlations between different detectors (off diagonal elements) decrease at higher frequencies. In the case of polarbear-i, polarization 
noise turned out to be fairly white in the chosen science bands, for 500 < i < 2000. The squared patterns noticeable in the left panel 
correspond to common structures of the focal plane: detectors are physically close and share some parts of the readout, leading to expected 
correlation. 


construction of , based on the stationarity assumption. 

In addition, similarly to the computation of in Eq. (25), accelerating or decelerating motions of the telescope 
(turnarounds) are also i nclu ded, as they did not show to have significantly different noise properties. We use the 
symmetry relation, Eq. (B4), to speed up the computation time. 

Besides, bins are chosen to be logarithmically spaced in the [0.01, 15] Hz range. An interesting possibility is to 
use an adaptive binning, which would make sure that any important feature in the noise power spectra is optimally 
encapsulated by the covariance matrix. 


TOD simulation given 

In this paragraph we propose a procedure to simulate realistic TODs using the informations encapsulated in the 
quantity. We follow the usual method to simulate noise-only TOD from the full binned noise covariance matrix, D. 
Eirst, we define 

V b, L,5 ^ ybj (B5) 

where is any operation such that, for any A which has a square root, 

VaVa^ = A (B6) 

Eor each frequency /, we generate a ridet random vector G We can write 

+ (B7) 

where and G M are independent random variables. The simulated noise time stream, n/, expressed in Eourier 
domain, will be given by 


V/ e {/„,in,/max}, n/ = ^ 


(B8) 


in which equation one should remember that b = b{f). In addition, /min = 1/^, T being the length of the simulated 
TOD, and /max is the maximum frequency achievable. In our implementation we choose T = (length of a CES), 
corresponding to /min < min(6). Therefore, we decide to apply a same power to the lowest frequency modes of 

the TODs. Of course, one would like ha ve a first bin b as low as possible in frequency. However, one has to take into 
account the average over S used in Eq. (B3): taking S ^ 2b consecutive subscans, leading to ^ 5 — 10 groups of S in 
the case of polarbear-i scanning strategy, turned out to give enough informations about the lowest frequency modes 
of the TODs. 

It is possible to interpolate L-^- over the frequency bins, leading to a L-^ matrix. Einally, the time domain simulated 
noise, for a detector i, is therefore given by 


= EET 


(»/) . 


m 
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In practice, to implement Eq. (B5), we diagonalize D/- 


ij' 


V b, = E ® Vk^ 

k 


(BIO) 


where and Vk are respectively the eigen values and eigen vectors. We then simply write 


V b, Lj5 = y] \/^ Vk ® Vk^ 
k 


(Bll) 


Besides, our implementation uses a simple parallelization scheme, each pro cess doing 


CES, i.e., Eqs. (B8) and (B9). All the 
for the entire season. 


and their square root (Eq. (B5)), are 


the TODs-simulation of a given 
precomputed and saved on disk 


Discussion 

Eor a given frequency bin 6, the matrix has a rank R = x where is the number of frequencies 
included in the bin b and is the number of groups of subscans. It turns out that, for a given 6, is not full rank 
i.e., R < Consequently, there are — R singular eigen values, making D-^ semi-positive definite for a given 

b —and preventing for example its Cholesky decomposition. We checked that this singularity did not have any impact 
on the simulations. 
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